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1  -  INTRODUCTION 

This  report  summarizes  the  results  of  the  SERDP  SEED  Project  UX-1284,  Application  of 
Wavelets  for  Detection  and  Discrimination  of  UXO.  This  project  addressed  the 
Statement  of  Need  (SON)  for  SERDP’s  SEEDSON-02-04,  which  calls  for  improved 
signal  processing  and  or  sensors  to  aid  in  discrimination  of  clutter  from  ordnance  in 
contaminated  areas  where  overlapping  signatures  are  common.  We  investigated  the 
effectiveness  of  signal  processing  techniques  based  on  wavelets  (1)  to  improve  the  signal 
to  noise  ratio  and  extract  additional  information  from  the  signals  and  (2)  as  part  of  a 
probability-based  approach  for  discriminating  between  ordnance  and  clutter. 

Previous  work  with  wavelet-based  signal  processing  has  shown  the  ability  of  wavelet- 
based  filtering  techniques  to  significantly  increase  the  signal  to  noise  ratio,  allowing 
weak,  but  real,  signals  to  be  extracted  from  noisy  backgrounds.  The  increase  in  the  signal 
to  noise  ratio  is  also  expected  to  result  in  a  more  accurate  and  detailed  description  of 
magnetic  anomalies  that  are  currently  detectable.  This  increase  in  available  detail  may 
make  it  easier  to  distinguish  between  multiple  overlapping  anomalies. 

Improved  methods  for  detection  and  discrimination  of  UXO  are  desirable  because  the 
size  of  the  ordnance  contaminated  areas  are  large  -  estimated  by  the  U.S.  DoD  to  be 
approximately  1 1  million  acres  (roughly  equal  to  the  combined  areas  of  the  states  of 
Vermont  and  New  Hampshire).  If  successful,  this  demonstration  will  make  the  accuracy 
of  reconnaissance  of  unexploded  ordnance  using  airborne  methods  more  comparable  to 
ground-based  detection  methods,  and  provide  a  method  for  improving  the  detection  and 
discrimination  capabilities  of  ground-based  methods. 

This  report  is  organized  as  follows:  Section  2  presents  background  information  on 
wavelets,  wavelet  decomposition  using  dyadic  wavelets  and  wavelet  packets,  and 
describes  the  application  of  probabilistic  methods  for  classification.  Section  3  describes 
the  project  activities  associated  with  applying  wavelet  filtering  to  airborne  magnetic  data. 
The  selection  of  the  wavelet-filtering  scheme  is  described  and  results  from  wavelet 
filtering  are  presented  and  compared  to  standard  filtering  results.  Section  4  describes  the 
project  activities  associated  with  discrimination  of  ordnance  from  clutter  using 
probabilistic  methods.  The  main  focus  of  this  effort  was  limited  to  wavelet-based 
descriptor  selection.  Section  5  presents  conclusions  and  recommendations  based  on  the 
project  results. 


Report  Documentation  Page 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

14  AUG  2003 


2.  REPORT  TYPE 


3.  DATES  COVERED 

00-00-2003  to  00-00-2003 


4.  TITLE  AND  SUBTITLE 

SERDP  SEED  Project  UX-1284:  Application  of  Wavelets  for  Detection 
and  Discrimination  of  UXO 

6.  AUTHOR(S) 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Oak  Ridge  National  Laboratory, Oak  Ridge, TN, 37831-6010 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


5a.  CONTRACT  NUMBER 


5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

5d.  PROIECT  NUMBER 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 


13.  SUPPLEMENTARY  NOTES 


14.  ABSTRACT 


15.  SUBIECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 


a.  REPORT 

unclassified 


b.  ABSTRACT 

unclassified 


c.  THIS  PAGE 

unclassified 


17.  LIMITATION  OF 

18.  NUMBER 

ABSTRACT 

OF  PAGES 

Same  as 

25 

Report  (SAR) 

19a.  NAME  OF 
RESPONSIBLE  PERSON 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


2  -  BACKGROUND  INFORMATION 


2.1  Wavelet  Transform  Description 

In  wavelet  decomposition,  the  signal  is  expanded  in  a  set  of  orthonormal,  compactly 
supported  basis  functions.  This  property  allows  wavelets  to  preserve  both  time  (spatial) 
and  frequency  information.  The  ability  to  preserve  time  information  makes  the  wavelet 
transform  particularly  useful  for  describing  transient  or  non-stationary  data  and  for 
describing  data  sets  containing  discontinuities. 


The  orthonormal  basis  used  in  the  wavelet  decomposition  is  based  on  a  single  function 
yKt)  that  is  defined  recursively  by  dilation  and  scaling 


Va.b 


t  -  b  N 
a  ; 


(1) 


where  a  specifies  the  wavelet  scale  and  b  specifies  the  wavelet  translation.1  The 
continuous  wavelet  transform  of  a  signal  x(t)  is  defined  as 

xv(a,b)  =  ]V0/ O  x(t)dt.  (2) 

~oo 

The  continuous  wavelet  transform  can  be  extended  to  a  discretely  sampled  data  set  x(t)  of 
N  samples  defined  as 


x(t)  =  x„  =  x(tn),  tn  =  nAt,  n  =  1,2 . N ,  (3) 

where  At  is  the  sampling  interval.  It  can  be  shown  that  if  the  scale  a  and  the  shift  b  are 
discretized  by  using 

a  -  2'j  At  ,  b  =  —  At ,  (4) 

2' 


where  j  and  k  are  integers,  an  orthonormal  basis  function  y/jk  can  be  derived; 

l 

¥j* =  ^hy/(  2in  ~ k)'  (5) 

The  discrete  wavelet  transform  is  then  given  by  the  summation 

djik  =  Va7  X  2L2W  (  2j  n  -  k  )  xn,  (6) 

n~0 

where  the  wavelet  transform  coefficients  are  denoted  by  djt k-  The  wavelet  function  is 
defined  by  a  set  of  quadrature  mirror  filter  coefficients.  Furthermore,  Daubechies  has 
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shown  that  these  filter  coefficients  can  be  selected  such  that  the  resulting  wavelet  has  M 
vanishing  moments,  that  is 


J  y/(  x  )  xm  d  x  =  0,  m  =  0, 1 . .  M  -  1 . 2  (7) 

Wavelet  functions  with  compact  support  and  vanishing  moments  result  in  local  events 
affecting  only  a  limited  number  of  wavelet  coefficients  and  in  smooth  portions  of  the 
signal  being  accurately  represented  by  few  coefficients.  These  properties  are  useful  for 
data  compression  and  for  denoising  signals  because  virtually  all  the  significant 
information  is  contained  in  relatively  few  wavelet  coefficients.  The  remaining 
coefficients  contain  noise  and  can  be  set  equal  to  zero,  effectively  eliminating  the  noise. 


The  application  of  the  wavelet  decomposition  (or  wavelet  filter)  results  in  two  subbands, 
each  with  half  the  bandwidth  of  the  original  signal.  The  high-frequency  subband  is 
generally  referred  to  as  "signal  details";  if  the  highest  frequency  contained  in  the  original 
signal  is  a>hi,  the  frequency  range  in  this  subband  extends  from  0)hi/ 2  to  0)hj .  The  low- 


free  Decomposition 


frequency  subband  is  referred  to  as  the 
"signal  approximation"  and  has  a 
frequency  range  of  0  to  cohJ 2 .  Each 
subband  has  half  the  number  of  samples  as 
the  original  signal. 

In  the  dyadic  wavelet  transform,  the 
wavelet  decomposition  is  recursively 
applied  to  the  signal  approximation  a 
specified  number  of  times.  Each 
application  of  the  wavelet  decomposition 
results  in  formation  of  a  wavelet  level, 
with  each  level  composed  of  a  detail 
subband  and  an  approximation  subband. 
The  approximation  subband  is  then 
operated  on  again  by  the  wavelet 
decomposition,  forming  the  next  wavelet 
level. 

Figure  1  shows  a  3-level  dyadic  wavelet 
transform  tree.  If  the  original  signal  has  N 
samples,  the  first  level  detail  subband, 
identified  as  <71 ,  contains  N/2  elements, 
the  second  level  detail  subband,  <72, 
contains  N/ 4  elements,  and  the  third  level 
each  contain 

samples,  <71  has  a 


detail  subband,  <73,  and  third  level  approximation  subband,  a3, 

N/ 8  elements.  If  the  original  signal  has  time  increment  of  At  between 
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time  increment  of  2At  between  samples,  dl  has  a  time  increment  of  4 At  between  samples, 
and  d3  and  a3  have  a  time  increment  of  8  At  between  samples.  In  an  m-element  dyadic 
wavelet  transform,  the  wavelet  decomposition  is  applied  m  times  and  each  detail  subband 
contains  half  as  many  elements  as  the  previous  detail  subband.  Table  1  shows  the 
number  of  elements,  the  frequency  range,  and  the  time  between  elements  for  each  of  the 
subbands  in  the  3-level  dyadic  wavelet  transform. 


2.2  Wavelet  Packet  Description 


subband 

number  of 
elements 

frequency  range 

time  increment 
between  elements 

original  signal 

N 

0-  CQhi 

At 

dl 

N/2 

“>Hi/2  -  <°Hi 

2  At 

dl 

N/4 

*W4  -  03 Hi/ 2 

4  At 

d3 

N/  8 

ffl/,,/8  -  0)hi/4 

8At 

a3 

N!  8 

0-  0)hi/8 

OO 

Table  1.  Number  of  elements,  frequency  range,  and  time  increment  for  the 
subbands  of  the  3-element  dyadic  wavelet  transform. 


The  wavelet  packet  is  much  like  the  dyadic  wavelet  transform,  except  the  wavelet 
decomposition  is  applied  to  the  detail  subbands  as  well  as  to  the  approximation  subbands. 
Figure  2  shows  a  3-level  wavelet  packet  tree.  The  result  of  the  wavelet  packet 
decomposition  is  the  subbands  a3,  cf3\,  d32,  d3 3,  d3 4,  d35,  d36,  and  d31 .  In  an  m- 
element  wavelet  packet  transform,  the  wavelet  decomposition  is  applied  2m  times  and  all 
subbands  contain  the  same  number  of  elements,  have  the  same  time  increment  between 
elements,  and  have  equal  bandwidths.  Table  2  shows  the  number  of  elements,  the 
frequency  range,  and  the  time  increment  between  elements  for  each  wavelet  packet 
subband. 

Comparison  of  the  dyadic  wavelet  transform  and  the  wavelet  packet  transform  shows  that 
the  wavelet  packet  transform  requires  considerably  more  computations  and  results  in  a 
finer  frequency  decomposition  and  a  coarser  time  decomposition. 

The  choice  of  using  either  the  dyadic  wavelet  transform  or  the  wavelet  packet  transform 
will  obviously  be  problem-dependent.  Wavelet  filtering  of  magnetometer  data  requires  a 
relatively  fine  frequency  decomposition  in  order  to  remove  helicopter  rotor  noise  without 
seriously  affecting  the  frequency  range  containing  UXO  anomaly  information.  For  this 
reason,  we  have  selected  the  wavelet  packet  transform  for  wavelet  filtering. 
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Figure  2.  Three  level  wavelet  packet  tree. 


subband 

number  of 
elements 

frequency  range 

time  increment 
between  elements 

original  signal 

N 

0  - 

At 

<m 

N/  8 

®*,/8  -  0Jhl/4 

SAt 

d32 

Nl  8 

* W 4  -  3fflw/8 

SAt 

d33 

Nl  8 

3®„,/8  -  CO.J2 

SAt 

d34 

Nil 8 

®«/2  -  5(Ohi/S 

SAt 

d3  5 

NIS 

5ahi/S  -  3(0 hi /4 

SAt 

d36 

Nf  8 

3©w/4  -  7®w/8 

SAt 

d31 

NIS 

7®w/8  - 

8  At 

a\ 

JW8 

0  -  ©w/8 

SAt 

Table  2.  Number  of  elements,  frequency  range,  and  time  increment  for  the  subbands 
of  the  three-element  wavelet  packet  transform. 
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2.3  Application  of  Probabilistic  Methods  for  Classification 

Applying  probabilistic  methods  for  classification  involves  several  steps.  First,  a  set  of 
descriptors  that  adequately  distinguish  between  classes  must  be  identified.  Next,  the 
statistical  properties  of  the  descriptor  set  for  each  class  must  be  determined,  usually  by 
measuring  the  descriptors  for  samples  taken  from  members  of  each  class  and  calculating 
the  usual  statistical  quantities  such  as  the  average  and  standard  deviation.  Finally,  a 
statistical  test  is  developed  using  the  descriptor  statistical  properties  to  determine  class 
membership  of  an  unknown  sample. 

Our  original  intention  in  this  project  was  to  perform  all  three  steps  for  a  limited  number 
of  samples  of  both  ordnance  and  clutter.  Unfortunately,  because  of  unexpected  difficulty 
in  the  performance  of  the  wavelet  filtering  portion  of  the  project,  we  were  forced  to  focus 
only  on  the  identification  of  the  descriptor  set.  It  has  been  our  experience  that  if  an 
adequate  descriptor  set  can  be  identified,  that  is,  one  that  can  distinguish  between  the 
classes,  then  the  other  two  steps  involved  in  applying  probabilistic  classification  methods 
will  present  little  difficulty.  By  focusing  our  effort  on  the  key  step  in  the  process,  we  feel 
we  can  most  effectively  evaluate  the  potential  for  applying  probabilistic  methods  for 
classification  of  ordnance  and  clutter,  given  the  project  constraints. 


3  -  WAVELET  FILTERING  OF  AIRBORNE  MAGNETOMETER  DATA 
3.1  Approach 

The  wavelet  filtering  should  retain  signal  features  containing  magnetic  anomaly 
information  while  removing  noise  and  signal  features  describing  other  phenomena  such 
as  rotor  noise,  large-scale  magnetic  effects,  etc.  The  approach  used  to  determine  the 
optimum  wavelet-filtering  scheme  was  as  follows: 

1)  Generate  synthetic  data  for  several  magnetic  anomalies  using  the  MAGMOD 
computer  code.  Three  different  ordnance  types  buried  at  three  different  depths  were 
used,  resulting  in  9  synthetic  magnetic  anomalies. 

2)  Extract  a  pure  noise  signal  from  the  measured  magnetometer  data.  A  portion  of  the 
magnetometer  signal  collected  at  an  altitude  greater  than  1 0  meters  was  used  as  the  noise 
signal.  Experience  shows  that  magnetic  anomalies  indicative  of  UXO  are  not  detectable 
at  altitudes  greater  than  10  meters,  thus  the  magnetometer  measurements  at  altitudes 
greater  than  10  meters  contains  only  effects  such  as  rotor  noise  and  large-scale  magnetic 
anomalies. 

3)  Apply  the  wavelet  packet  decomposition  to  each  signal  and  calculate  the  energy 
contained  in  each  subband.  Ideally,  we  would  find  no  overlap  in  subband  signal  energy 
between  the  two  signals. 
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4)  The  optimum  wavelet-filtering  scheme  would  retain  the  subbands  containing  high 
synthetic  signal  energy  and  discard  subbands  containing  high  noise  signal  energy.  The 
filtered  signal  is  formed  by  inverse  wavelet  transforming  the  retained  subbands. 

5)  The  effectiveness  of  the  wavelet-filtering  scheme  can  be  evaluated  by  comparing  the 
wavelet  filtering  results  to  standard  filtering  results.  Although  wavelet  filtering  will  be 
performed  using  profile  data,  the  comparison  can  also  be  performed  after  converting  the 
signals  to  gridded  data.  The  wavelet  filtering  would  be  considered  superior  to  standard 
filtering  if 

a)  additional  anomaly  detail  can  be  discerned  in  the  wavelet  filtered  results, 

b)  additional  anomalies  can  be  detected  in  the  wavelet  filtered  results, 

c)  the  overall  noise  level  is  reduced  in  the  wavelet  filtered  results,  or 

d)  the  application  of  wavelet  filtering  is  significantly  faster  or  requires 

significantly  less  user  intervention  during  the  filtering  process. 

The  wavelet  analysis  is  described  in  the  next  section.  Descriptions  of  the  synthetic  and 
noise  signals  are  given  and  the  results  of  the  wavelet  packet  analysis  for  each  signal  are 
presented.  The  selection  of  the  wavelet-filtering  scheme  is  described  and  the  results  of 
this  scheme  applied  to  the  synthetic  signal  with  added  noise  are  shown.  Section  3 
describes  the  magnetometer  data  used  in  to  compare  the  filtering  and  shows  the  results  of 
the  wavelet-filtering  scheme  on  the  collected  data.  A  comparison  of  the  wavelet-filtered 
and  standard-filtered  data  is  given  for  both  profile  and  gridded  data.  Conclusions  and  a 
summary  of  the  analysis  results  are  contained  in  the  final  section. 


3.2  Wavelet  Filtering  Scheme  Selection 

The  wavelet-filtering  scheme  was  selected  after  examining  the  energy  in  the  wavelet 
subbands  for  the  synthetic  signal,  which  contains  only  anomaly  signatures,  and  a  pure 
noise  signal.  For  each  signal,  a  four  level  wavelet  packet  decomposition  was  performed 
on  the  synthetic  data  using  a  Daubechies-20  coefficient  wavelet. 

3.2.1  Synthetic  Signal 

The  MAGMOD  computer  code  was  used  to  calculate  the  magnetometer  responses  of  for 
three  different  ordnance  types  buried  at  three  different  depths.  The  synthetic  signal  is 
shown  in  Figure  3. 

The  sample  number  corresponds  to  the  sensor  traveling  from  south  to  north;  that  is,  the 
first  data  point  represents  the  southernmost  point.  There  are  9  different  magnetic 
anomalies  from  bodies  centered  every  50  m  from  x  =  0,  i.e.  at  x  =  (0,50, 100, ...400).  The 
he  sample  number  corresponds  to  the  sensor  traveling  from  south  to  north;  that  is,  the 
first  data  point  represents  the  southernmost  point.  There  are  9  different  magnetic 
anomalies  from  bodies  centered  every  50  m  from  x  =  0,  i.e.  at  x  =  (0,50, 100,. ..400).  The 
three  southernmost  anomalies  are  from  76-mm  ordnance  1.5, 2.5,  and  3.5  m  from  the 
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Figure  3.  The  MAGMOD-calculated  synthetic  magnetic  anomaly  signal. 

sensor.  The  next  three  anomalies  are  from  105-mm  ordnance  1.5, 2.5,  and  3.5  m  from  the 
sensor.  The  three  northernmost  anomalies  are  from  155-mm  ordnance  1.5, 2.5,  and  3.5  m 
from  the  sensor.  These  depths  correspond  to  a  flight  height  of  1 .5  m  above  ground  level 
and  burial  depths  of  0,  1,  and  2  m.  Each  object  is  angled  45  degrees  below  horizontal, 
pointed  north-south  with  the  southernmost  end  of  the  ordnance  nearest  the  surface.  The 
ambient  magnetic  field  is  56000  nT,  has  a  declination  of  0  degrees,  and  an  inclination  of 
70  degrees.  The  helicopter  velocity  was  assumed  to  be  10  m/sec  and  the  sample  rate  was 
60  Hz. 

The  four-level  wavelet  packet  tree  representing  the  wavelet  decomposition  of  the  signal  is 
shown  in  Figure  4.  Each  branch  point  represents  a  time  series  obtained  from  wavelet 
decomposition.  For  example,  the  original  time  series,  represented  by  (0,0),  is 
decomposed  into  two  subbands,  (1,0)  and  (1,1).  The  (1,0)  time  series  is  then  decomposed 
into  subbands  (2,0)  and  (2,1).  This  process  is  repeated  to  form  the  wavelet  packet  tree 
shown  in  Figure  4.  The  time  series  at  the  end  of  each  "branch"  are  the  final  subbands 
used  to  represent  the  signal.  For  this  simple  tree,  the  original  signal  is  broken  up  into  the 
five  subbands  (1,1),  (2,1),  (3,0),  (4,2),  and  (4,3). 

After  the  signal  is  decomposed,  the  energy  in  each  subband  is  calculated.  Table  3  shows 
the  approximate  frequency  range  and  the  signal  energy  in  each  wavelet.  The  results  of 
the  energy  analysis  show  that,  besides  the  low  frequency  energy,  the  majority  of  the 
signal  energy  is  contained  in  subband  (4,2).  The  majority  of  the  signal  energy  would  be 
retained  if  subbands  (3,0)  and  (4,2)  are  used  to  reconstruct  the  signal. 
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Subband 


Approximate  frequency 
_ range  (Hz) _ 


Subband  Ei 
(%  Total  En 


(Lll 

i m 

(3,0) 

(4.2) 

(4.3) 


15-30 
7.5-15 
0-3.75 
3.75-5.625 
5.625  -  7.5 


0.003 

0.014 

95.25 

4.53 

0.22 


Table  3  .  Wavelet  subband  energy  for  synthetic  anomaly  signal. 
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3.2.2  Noise  Signal 


The  noise  signal  was  obtained  from  measured  data  taken  at  altitudes  above  10  meters.  At 
these  altitudes,  ground  magnetic  anomalies  have  a  minimal  influence  on  the  magnetic 
measurements;  the  sensors  respond  primarily  to  the  magnetic  field  of  the  helicopter  and 
to  fluctuations  in  the  earth's  magnetic  field.  The  pure  noise  signal  is  shown  in  Figure  5 
The  approximate  frequency  range  and  the  signal  energy  in  each  of  the  wavelet  subbands 


Noise  Signal 


5.2815 


5.2785 


5:278 


3  4 

sample  number 


Figure  5.  The  pure  noise  signal. 

are  shown  in  Table  4.  The  results  of  the  energy  analysis  show  that,  besides  the  low 
frequency  energy,  the  majority  of  the  signal  energy  is  contained  in  subbands  (1,1),  (2,1), 
and  (4,3).  If  these  subbands  are  not  included  in  the  reconstructed  signal,  the  majority  of 
the  noise  at  frequencies  greater  than  3.75  Hz  would  eliminated.  A  standard  regional 
correction  or  possibly  high  pass  filtering  could  remove  the  low  frequency  noise. 


Subband 

Approximate  frequency 
range  (Hz) 

M5ESM 

15-30 

2.66 

mkssm 

7.5-15 

1.27 

■ 

0  -  3.75 

94.22 

3.75  -5.625 

0.41 

■w 

5.625  -  7.5 

1.45 

Table  4  .  Wavelet  subband  energy  for  measured  noise  signal. 
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More  can  be  said  concerning  neglecting  the  low  frequency  components  to  eliminate  low 
frequency  noise.  Originally,  we  had  tried  using  a  longer  wavelet  decomposition  that 
broke  subband  (3,0)  and  additional  four  levels.  The  lowest  frequency  subband  contained 
a  significant  amount  of  noise  energy  and  was  neglected  during  signal  reconstruction.  The 
unexpected  result  was  the  appearance  of  a  relatively  large  scale,  low  frequency  artifact 
that  distorted  the  data  to  the  point  of  making  it  useless.  For  this  reason,  the  low 
frequency  subbands  have  been  retained  and  this  portion  of  the  standard  filtering  (i.e.,  the 
regional  correction)  was  used  to  remove  the  low  frequency  noise.  It  should  be  possible  to 
develop  a  wavelet-based  filtering  scheme  that  will  remove  the  low  frequency  noise 
without  adding  unacceptable  artifacts;  this  investigation  will  be  need  to  wait  until  the 
next  phase  of  the  follow-on  project. 

Additional  insight  into  the  performance  of  the  wavelet  decomposition  can  be  obtained 
from  the  information  shown  in  Figure  6.  In  this  figure,  the  magnitudes  of  the  discrete 
Fourier  transform  of  each  wavelet  subband  are  shown.  In  addition  to  the  obvious 
frequency  range  of  each  subband,  the  results  in  the  frequency  domain  show  the  rotor 
noise  at  approximately  6.5  Hz  and  an  unidentified  frequency  peak  at  5  Hz. 

Based  on  the  wavelet  subband  energy  in  the  synthetic  anomaly  signal  and  the  pure  noise 
signal,  the  wavelet  filtering  scheme  selected  is  to  use  only  subbands  (3,0)  and  (4,2)  in  the 


10  15  /Ll  .  20 

Frequency  (Hz) 


subband  (1,1) 
subband  (2,1) 
subband  (3,0) 
subband  (4,3) 
subband  (4,2) 


Figure  6.  Discrete  Fourier  transform  magnitude  of  the  wavelet  subbands  for  the 
pure  noise  signal. 
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signal  reconstruction.  Thus,  the  wavelet  filtering  will  consist  of  decomposing  the 
magnetometer  signal  by  using  the  wavelet  tree  structure  shown  in  Figure  4  (Daubechies 
20  coefficient  wavelet  is  used),  and  then  reconstructing  the  filtered  time  series  from 
subbands  (3,0)  and  (4,2).  This  filtering  scheme  passes  signals  in  the  frequency  range  of  0 
to  5.625  Hz,  effectively  eliminating  the  rotor  noise. 


3.3  Filtering  Results 

The  wavelet  filtering  strategy  was  applied  to  data  collected  from  the  Badlands  bombing 
range  in  September  2000.  The  data  was  collected  by  using  the  Oak  Ridge  Airborne 
Geophysical  System  -  Hammerhead  Array  (ORAGS™-HA).  Figure  7  shows  the  survey 
helicopter  with  the  mounted  ORAGS™-HA  system. 


Figure  7.  The  ORAGSTM  Hammerhead  array  for  UXO  detection. 

This  system  deploys  8  magnetometers  at  1.75m  with  three  magnetometers  located  in  each 
of  the  lateral  booms  and  two  magnetometers  located  in  the  forward  boom.  The  system 
uses  a  PC-based  data  acquisition  system  that  is  capable  of  recording  raw  synchronous 
Larmor  frequency  measurements  from  the  magnetometers  at  1200  Hz  sample  rate.  The 
ORAGS  system  can  record  raw  data  at  1200  Hz  sample  rate;  For  these  tests,  data  were 
downsampled  to  60Hz  in  order  to  minimize  the  impact  of  power  line  interference. 

Over  2200  acres  of  the  BBR  were  surveyed  with  the  ORAGS-HA  system  during  a  10  day 
period  in  September  2000.  This  included  a  test  site,  five  large  target  areas,  and  numerous 
transects.  The  test  site  and  target  areas  were  surveyed  as  part  of  an  Environmental 
Science  and  Technology  Certification  Program  (ESTCP)  demonstration  project,  while  the 
transects  were  surveyed  to  support  ongoing  Engineering  Evaluation/Cost  Analysis 
(EE/CA)  efforts  by  Parsons  Engineering  Science.  The  test  grid  was  100m  x  150m  and 
included  several  pieces  of  ordnance  that  are  not  common  at  BBR,  but  frequently  occur  at 
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other  UXO  sites.  These  include  105mm  and  155mm  artillery,  20mm  rounds,  mortars  and 
rockets.  Figure  8  shows  a  portion  of  a  typical  magnetometer  signal,  in  this  case  channel 
0,  as  the  helicopter  passes  over  a  magnetic  anomaly.  The  presence  of  two  magnetic 
anomalies  is  clearly  shown,  but  also  a  considerable  amount  of  noise  (primarily  rotor 
noise)  can  also  be  plainly  seen.  Three  profiles  are  shown.  The  first  shows  the  unfiltered 
data  sampled  at  1200  Hz.  These  data  have  undergone  the  compensation  correction,  but 
have  had  no  additional  signal  processing.  The  second  profile  shows  the  data  after 
standard  filtering.  In  this  case,  the  data  was  downsampled  to  60  Hz  before  filtering.  The 
results  show  a  significant  reduction  in  the  noise  level.  The  third  profile  shows  the  data 
after  wavelet  filtering.  Here,  the  filtered  data  was  upsampled  to  the  original  1200  Hz 
sampling. 


(a)  -  Unfiltered  signal,  sampled  at  1200  Hz. 


(b)  -  Results  of  current  filter  scheme  signal,  displayed  at  60  Hz 


(c)  -  Wavelet-filtered  signal,  displayed  at  1200  Hz  sample 
rate. 


Figure  8.  Conroarison  of  filter  results  on  magnetic  profile  data. 
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A  qualitative  comparison  of  the  filtered  results  shows  that  that  both  accurately  preserve 
the  shape  and  amplitude  of  the  peaks  associated  with  the  magnetic  anomaly.  However, 
the  wavelet-filtered  data  has  significantly  less  noise  than  the  standard-filtered  data.  Thus, 
it  appears  that  the  wavelet  filter  is  superior  to  the  standard  filter  for  removing  noise. 

Investigation  of  the  filtered  results  in  the  frequency  domain  provides  additional  insight 
into  the  performance  of  the  two  filters.  Figure  9  shows  the  filtered  results  in  the 
frequency  domain.  The  unfiltered  signal  shows  the  signal  energy  to  be  primarily  at  low 
frequencies  with  a  prominent  rotor  noise  peak  at  approximately  6.5  Hz  and  an 
unidentified  noise  peak  at  approximately  5  Hz.  Some  additional  signal  energy  is  present 
in  the  13  to  15  Hz  frequency  range.  After  standard  processing,  the  signal  energy  below 
approximately  4  Hz  is  preserved;  the  signal  content  above  4  Hz  has  been  completely 
eliminated.  The  standard  processing  effectively  removes  the  rotor  noise,  but  also 
removes  any  information  describing  the  magnetic  anomalies  between  4  and  6  Hz.  The 
unidentified  noise  peak  at  5  Hz  is  removed  by  the  standard  processing. 

After  wavelet-filtering,  the  signal  energy  is  preserved  up  to  a  frequency  of  approximately 
6  Hz.  The  rotor  noise  is  effectively  removed  and  additional  information  in  the  4  to  6  Hz 
frequency  range  is  retained.  This  additional  information  should  result  in  a  more  detailed 
representation  of  magnetic  anomalies  in  the  wavelet-filtered  data,  compared  to  the 
standard-filtered  data.  The  unidentified  noise  peak  at  5  Hz  is  preserved  in  the  wavelet- 
filtered  data. 

The  filtered  results  can  also  be  compared  after  gridding.  The  compensated  data  was 
filtered  using  both  the  standard  filtering  and  wavelet-filtering.  The  two  results  were 
gridded  by  using  GeoSoft  software.  The  results  are  shown  in  Figures  10  and  1 1.  Figure 
10  shows  a  portion  of  the  collected  data  for  a  relatively  large  area.  Close  examination  of 
the  figure  shows  the  wavelet-filtered  results  to  represent  the  magnetic  anomalies  with 
additional  detail  compared  to  the  standard-filtered  results.  The  region  enclosed  in  the 
circle  is  expanded  in  Figure  11.  Here,  the  additional  detail  preserved  by  the  wavelet 
filtering  is  clearly  evident.  Not  only  are  more  anomalies  shown,  but  the  anomalies  shown 
when  using  standard  filtering  are  shown  with  additional  detail  when  wavelet  filtering  is 
used.  This  additional  detail  can  be  attributed  to  the  additional  information  in  the  4  to  6 
Hz  frequency  range  that  is  preserved  by  wavelet  filtering. 
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FFT  Magnitude  FFT  Magnitude  FFT  Magnitude 


Raw  Signal  Frequency  Spectrum 


Figure  9.  Comparison  of  filtered  profile  data  in  the  frequency  domain. 
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Standard  filtered  Wavelet  filtered 


Figure  10.  Comparison  of  standard-filtered  and  wavelet  filtered  gridded 
data  -  laree  area  view. 


Standard  filtered  Wavelet  filtered 

Figure  11.  Comparison  of  standard-filtered  and  wavelet  filtered  gridded  data  - 
detailed  view. 
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A  quantitative  measure  of  the  difference  between  standard  filtering  and  wavelet  filtering 
can  be  obtained  from  analysis  of  high-altitude  noise  data.  The  measured  signal,  which  is 
collected  at  altitudes  greater  than  ten  meters,  can  be  filtered  using  each  filtering  method 
and  the  filtered  results  quantitatively  compared.  Figure  12  shows  the  frequency  domain 
representation  between  0  and  7  Hz  of  a  sample  of  high  altitude  data  after  filtering  using 
both  methods.  The  results  show  that  both  filtering  methods  are  practically  identical  at 
frequencies  below  4  Hz.  Between  4  and  6  Hz,  the  wavelet  filtering  preserves  more  of  the 
signal,  a  result  previously  shown  in  Figure  9.  Between  6  and  7  Hz,  the  frequency  band 
containing  rotor  noise,  it  is  seen  that  the  wavelet  filtering  removes  more  of  the  signal. 
Figure  13  shows  the  comparison  of  the  two  filtering  method  results  in  more  detail  over 
the  frequency  range  of  4  to  7  Hz. 


Figure  12.  Frequency  domain  comparison  of  high  altitude  noise  filtered  using 
standard  and  wavelet  filtering  -  0  to  7  Hz. 

A  measure  of  the  relative  energy  in  each  filtered  signal  can  be  made  by  integrating  the 
filtered  results  over  frequency.  The  ratio  of  the  energies  in  each  signal  indicates  the 
relative  magnitude  of  the  energy  content  of  the  filtered  signals.  Table  5  shows  the  ratio  of 
the  signal  energy  of  the  two  filtering  methods  over  the  frequency  ranges  of  0  to  4  Hz,  4 
to  6  Hz,  and  6  to  7  Hz.  The  results  show  that  the  two  filtering  methods  have  nearly 
identical  energies  over  the  frequency  range  of  0  to  4  Hz.  Wavelet  filtering  retains 
approximately  42%  more  signal  energy  over  the  frequency  range  of  4  to  6  Hz.  Over  the 
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frequency  range  of  6  to  7Hz,  wavelet  filtering  removes  34%  more  noise  than  standard 
filtering.  Thus,  for  this  case  (which  we  believe  to  be  representative  of  noise  in  the 
ORAGS™-HA  system),  wavelet  filtering  retains  42%  more  energy  in  the  4  to  6  Hz 
frequency  range  while  rejecting  34%  more  noise  in  the  6  to  7  Hz  frequency  range. 


Figure  13.  Frequency  domain  comparison  of  high  altitude  noise  filtered  using 
standard  and  wavelet  filtering  -  4  to  7  Hz. 


Frequency  Range 
(Hz) 

Signal  Energy  Ratio 
(wavelet  energy/standard  energy) 

%  Additional  Energy  in 
Wavelet-Filtered  Signal 

0-4 

1.004 

0.42% 

4-6 

1.425 

42.5% 

6-7 

0.66 

-34% 

Table  5.  Comparison  of  energy  content  of  the  filtered  signals. 
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A  comparison  of  the  ease  of  application  can  be  made  between  the  two  filtering  methods. 
Figure  14  compares  block  diagrams  of  the  filtering  steps  for  the  standard  filtering,  the 
current  wavelet  filtering,  and  our  expectation  of  the  future  wavelet-filtering  processes. 
Currently,  standard  filtering  requires  more  steps  then  the  wavelet  filtering.  Furthermore, 
some  of  the  filtering  steps,  such  as  comparing  regional  corrections  and  making  the  B- 
Spline  fit,  require  user  intervention  or  input.  This  intervention  may  be  relatively  benign 
when  filtering  small  data  sets,  but  can  significantly  increase  the  time  and  cost  for  filtering 
large  data  sets  and  such  intervention  requires  a  highly  skilled  analyst. 

Wavelet  filtering  greatly  simplifies  filtering  by  eliminating  several  of  the  current 
processing  steps,  including  some  that  require  user  intervention  or  input.  Future  wavelet 
filtering,  which  would  involve  converting  the  wavelet  filtering  into  a  GeoSoft  module, 


Current  Process 


Current  Wavelet  Process  Future  Wavelet  Process 


Figure  14.  Block  diagram  comparing  the  filtering  steps  for  standard  filtering, 
wavelet  filtering,  and  the  expected  final  wavelet  filtering  processes. 
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would  simplify  filtering  even  further,  possibly  making  it  essentially  a  "tum-key" 
operation  that  can  be  performed  without  user  interaction.  It  is  expected  that  the 
application  of  wavelet  filtering  can  significantly  decrease  the  manpower,  time,  and 
ultimately  the  cost  required  to  filter  large  data  sets.  Such  data  sets  will  become  common 
as  airborne  data  acquisition  is  used  to  survey  more  UXO  sites. 

3.4  -  Summary  of  Wavelet  Filtering  Results 

The  results  of  comparing  standard  filtering  with  wavelet  filtering  indicate  the  following: 

1)  Wavelet  filtering  removes  more  noise  at  frequencies  above  6  Hz  than  standard 
filtering.  This  result  is  qualitatively  shown  in  Figures  8, 12,  and  13.  For  the  analyzed 
noise  sample  it  was  found  that  wavelet  filtering  removed  approximately  34%  more  noise 
in  the  6  to  7  Hz  frequency  range  than  did  standard  filtering. 

2)  Wavelet  filtering  preserves  more  useful  information  in  the  4  to  6  Hz  frequency  range 
than  standard  filtering.  This  result  is  shown  in  the  comparison  of  the  frequency  content 
of  the  two  filtering  methods  shown  in  Figures  9,  12,  and  13  and  in  the  comparison  of  the 
gridded  data  shown  in  Figures  10  and  11.  For  the  analyzed  noise  signal  it  was  found  that 
wavelet  filtering  preserved  approximately  42%  more  signal  in  the  4  to  6Hz  frequency 
range  than  did  standard  filtering. 

3)  Wavelet  filtering  is  simpler  and  requires  less  user  intervention  than  current  filtering.  It 
is  expected  that  wavelet  filtering  can  be  further  simplified  by  incorporating  the  filtering 
into  a  GeoSoft  module. 


4  -  DISCRIMINATION  OF  ORDNANCE  AND  CLUTTER 

As  described  above,  we  were  limited  in  this  project  to  identifying  a  partial  descriptor  set 
for  use  in  determining  if  a  magnetic  anomaly  represents  ordnance  or  clutter.  This 
descriptor  set  identification  was  done  with  an  extremely  limited  number  of  measured 
samples  of  ordnance  and  clutter.  The  samples  used  in  this  evaluation  were  two  practice 
bombs,  a  barbed  wire  fragment,  a  bomb  fin,  and  magnetically  responsive  soil.  Two 
channels  of  measured  data  were  used  in  the  evaluation  of  each  sample. 

4.1  Descriptor  Selection  and  Evaluation 

Our  efforts  focused  on  evaluating  the  potential  of  three  algorithmically-derived 
descriptors;  wavelet  subband  energy,  signal  entropy,  and  the  discrete  Fourier  transform 
spectra.  Wavelet  subband  energy  measures  the  amount  of  signal  energy  in  a  particular 
wavelet  subband.  This  quantity  indicates  the  frequency  distribution  of  energy  in  a  signal. 
It  is  reasonable  to  expect  that  the  energy  distribution  of  relatively  large,  intact  metallic 
bodies  such  as  ordnance  will  have  different  energy  distributions  than  smaller,  distributed 
metallic  bodies.  Signal  entropy  is  a  measure  of  correlation  or  noise  in  a  signal.  Again,  it 
is  expected  that  ordnance  signals  would  have  relatively  highly  correlated  signals, 


20 


resulting  in  relatively  low  entropy.  Clutter  signals  should  be  less  correlated,  resulting  in 
relatively  high  entropy.  The  discrete  Fourier  transform  is  another  measure  of  the 
distribution  of  signal  energy  over  frequency.  A  comparison  of  the  descriptor  values  for 
each  of  the  samples  is  presented  in  this  section. 

4.1.1  Wavelet  Subband  Energy 

A  four  level  dyadic  wavelet  decomposition  was  used  on  the  sample  signals.  An  8- 
coefficient  Daubechies  wavelet  was  used  in  the  decomposition.  The  wavelet  tree 
structure  describing  the  decomposition  is  shown  in  Figure  15.  The  energy  of  subband 
(4,1)  was  used  as  the  descriptor.  This  subband  corresponds  to  an  approximate  frequency 
range  of  1 .5  to  4  Hz.  Figure  16  shows  the  subband  energy  as  a  percentage  of  the  total 
energy  in  the  signal.  The  sample  number  corresponds  to  the  two  data  channels  for  the 
two  practice  bombs  and  the  two  data  channels  for  the  bomb  fin,  barbed  wire  sample,  and 
magnetically  responsive  soil.  The  results  show  that  in  general  the  practice  bombs  have  a 
much  higher  percentage  of  signal  energy  in  subband  (4,1)  than  the  bomb  fin,  barbed  wire 
sample,  or  the  magnetically  responsive  soil.  The  one  exception  is  one  barbed  wire 
sample  that  contains  an  energy  equivalent  to  the  practice  bomb. 


Figure  15.  Four-level  dyadic  wavelet  tree  structure  used  in 
the  wavelet  decomposition 


21 


Comparison  of  Signal  Energies  -  Wavelet  Subband  (4-1) 


2  3 

sample  number 


□  Barbed  Wire  ■  Hot  Dirt  m  Bomb  Fin  n  Bomb 


Figure  16.  Signal  energies  for  wavelet  subband  (4,1). 

4.1.2  Signal  Entropy 

The  signal  entropy  for  the  descriptor  samples  is  shown  in  Figure  17.  As  expected  the 
entropy  for  the  practice  bombs  are  considerably  less  than  that  of  the  clutter  objects. 
Again,  the  one  exception  is  one  of  the  barbed  wire  samples,  which  has  a  similar  entropy 
value  to  the  practice  bombs. 


Comparison  of  Signal  Entropies 


sample  number 


□  Barbed  Wire  ■  Hot  Dirt  H  Bomb  Fin  u  Bomb 


Figure  17.  Signal  entropies. 
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4.1.3  Discrete  Fourier  Transform 

The  discrete  Fourier  transform  is  another  measure  of  the  energy  distribution  of  a  signal. 
Figure  18  shows  the  Fourier  spectra  of  the  descriptor  samples.  To  generate  the  spectra 
shown  in  this  figure,  each  of  the  sample  signals  was  first  normalized  to  have  a  peak 
response  magnitude  of  ±  1.0.  is  shown.  Comparing  the  results,  one  sees  that  the  spectra 
of  the  ordnance  and  clutter  examples  are  quite  different.  In  general,  the  spectra  of  the 
ordnance  are  relatively  flat  with  no  prominent  peaks  in  the  0  to  4  Hz  frequency  range. 

The  spectra  of  the  clutter  samples  show  significant  peaks  in  the  0  to  4  Hz  frequency  range 
and  contain  less  energy  (as  a  percentage  of  total  energy)  in  the  2  to  4  Hz  frequency  range. 
Although  not  done  in  this  investigation,  it  should  be  possible  to  select  additional 
descriptors  from  the  discrete  Fourier  transform  to  distinguish  between  ordnance  and 
clutter. 


Figure  18.  Discrete  Fourier  transforms  spectra. 
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5  -CONCLUSIONS  AND  RECOMMENDATIONS 


Based  on  our  results,  we  believe  the  following  conclusions  are  justified: 

1)  Wavelet  filtering  removes  more  noise  at  frequencies  above  6  Hz  than  standard 
filtering.  This  result  is  qualitatively  shown  in  Figures  8, 12,  and  13.  For  the  analyzed 
noise  sample  it  was  found  that  wavelet  filtering  removed  approximately  34%  more  noise 
in  the  6  to  7  Hz  frequency  range  than  did  standard  filtering. 

2)  Wavelet  filtering  preserves  more  useful  information  in  the  4  to  6  Hz  frequency  range 
than  standard  filtering.  This  result  is  shown  in  the  comparison  of  the  frequency  content 
of  the  two  filtering  methods  shown  in  Figures  9,  12,  and  13  and  in  the  comparison  of  the 
gridded  data  shown  in  Figures  10  and  11.  For  the  analyzed  noise  signal  it  was  found  that 
wavelet  filtering  preserved  approximately  42%  more  signal  in  the  4  to  6Hz  frequency 
range  than  did  standard  filtering. 

3)  Wavelet  filtering  is  simpler  and  requires  less  user  intervention  than  current  filtering. 
Figure  12  compares  the  application  of  the  two  methods.  It  is  expected  that  wavelet 
filtering  can  be  further  simplified  by  incorporating  the  filtering  into  a  GeoSoft  module. 
This  reduction  in  complexity  would  result  in  significant  savings  in  both  time  and  cost, 
especially  if  large  areas  are  to  be  surveyed. 

4)  The  very  preliminary  work  done  during  this  study  encourages  us  that  a  descriptor  set 
that  can  be  used  to  distinguish  between  ordnance  and  clutter  can  be  identified.  This 
descriptor  set  would  most  likely  include  the  descriptors  used  in  this  study  but  would 
probably  include  additional  descriptors.  If  such  a  descriptor  set  can  be  identified,  it 
should  be  possible  to  use  these  descriptors  with  probabilistic  methods  to  classify  a 
magnetic  anomaly  as  being  either  ordnance  or  clutter. 

We  recognize  that  the  number  of  ordnance  and  clutter  samples  prevents  a  statistically 
valid  determination  on  the  ability  of  the  descriptor  set  to  distinguish  between  ordnance 
and  clutter.  However,  because  the  limited  results  are  so  unambiguous,  we  feel  confident 
in  concluding  that  there  is  real  potential  for  identifying  a  descriptor  set  that  will 
accurately  distinguish  between  ordnance  and  clutter. 

Based  on  our  results,  we  believe  that  both  wavelet  filtering  and  the  application  of 
probabilistic  methods  for  distinguishing  ordnance  from  clutter  have  potential  for 
improving  the  detection  and  identification  of  UXO.  We  recommend  that  a  follow-on 
project  be  initiated  to 

1)  Fully  develop  the  wavelet  filtering  technique  and  incorporate  the  final  wavelet  filtering 
into  a  GeoSoft  module.  This  module  would  be  made  available  to  the  UXO  community 
and  should  result  in  the  extraction  of  more  useful  information  from  the  measured  data 
than  is  currently  being  obtained.  In  addition  to  the  improvement  in  filtering  performance, 
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we  believe  that  the  use  of  wavelet  filtering  will  reduce  the  time,  manpower,  and  cost  of 
filtering  the  large  data  sets  typically  collected  by  airborne  UXO  surveys. 

2)  Fully  explore  the  potential  of  wavelet-based  descriptors  and  probabilistic  methods  for 
discriminating  between  ordnance  and  clutter.  The  preliminary  work  shown  here  appears 
promising  but  is  only  a  start  in  evaluating  the  potential  of  this  approach.  Evaluation 
using  larger  data  sets  involving  many  more  examples  of  ordnance  and  clutter  objects  is 
needed  to  determine  if  an  effective  descriptor  set  can  be  identified.  If  successful,  this 
approach  would  form  the  basis  for  a  screening  of  magnetic  anomalies,  allowing  analysts 
to  focus  on  only  the  most  promising  anomalies.  Again,  a  reduction  in  the  time, 
manpower,  and  cost  of  evaluating  the  survey  results  would  be  achieved. 


REFERENCES 

1 .  Akansu,  A.  N.,  and  R.  A.  Haddad,  Multiresolution  Signal  Decomposition,  1992,  Academic 
Press,  Inc.,  San  Diego,  CA.,  pp.  292-313. 

2.  Daubechies,  I.,  Ten  Lectures  on  Wavelets,  1992,  Society  for  Industrial  and  Applied 
Mathematics,  Philadelphia,  pp.  199-209. 


25 


